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Abstract: We propose a model for the propagation of a heavy-quark in a hot plasma, to be 
viewed as a first step towards a full description of the dynamics of heavy quark systems in a 
quark-gluon plasma, including bound state formation. The heavy quark is treated as a non 
relativistic particle interacting with a fluctuating field, whose correlator is determined by 
a hard thermal loop approximation. This approximation, which concerns only the medium 
in which the heavy quark propagates, is the only one that is made, and it can be improved. 
The dynamics of the heavy quark is given exactly by a quantum mechanical path integral 
that is calculated in this paper in the Euclidean space-time using numerical Monte Carlo 
techniques. The spectral function of the heavy quark in the medium is then reconstructed 
using a Maximum Entropy Method. The path integral is also evaluated exactly in the case 
where the mass of the heavy quark is infinite; one then recovers known results concerning 
the complex optical potential that controls the long time behavior of the heavy quark. The 
heavy quark correlator and its spectral function is also calculated semi-analytically at the 
one-loop order, which allows for a detailed description of the coupling between the heavy 
quark and the plasma collective modes. 

Keywords: Thermal Field Theory, Heavy Quark Physics. 
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1. Introduction 

Understanding the dynamics of heavy quarks in a quark-gluon plasma, and the fate of their 
possible bound states, has remained a difficult issue, ever since the original proposal of Mat- 
sui and Satz [1] to view the dissolution of J/^'s mesons in a quark-gluon plasma as a signal 
of deconfinement (for a recent review see for instance Ref. [2]). Aside from many studies 
based on the assumption that the dominant effect of the plasma is to screen the interaction 
potential, more recently, the problem has been attacked using a "first principle" approach, 
namely by calculating the QQ spectral functions reconstructed from the corresponding Eu- 
clidean correlators provided by lattice QCD. The melting of the J/V'j for instance, is then 
signaled by the disappearance of the corresponding peak in its spectral function. The first 
results of such an analysis led to the surprising result that the J /ip appears to survive till 
temperatures well above Tc [3, 4, 5, 6, 7], in sharp contrast with studies based on screened 
potentials. A comparison between correlators and spectral functions evaluated on the lat- 
tice and within different potential models was attempted in [8, 9, 10], revealing ambiguities 



in the whole procedure. Another hne of first principle calculations was undertaken in a 
number of recent papers [11, 12, 13, 14, 15, 16]: in these works, the correlator of a heavy 
quark pair is calculated directly in real time, revealing that the long time behavior of the 
dynamics can be encompassed by a Schrodinger equation with a complex potential that 
describes both the effects of screening and, through its imaginary part, of the collisions 
with the plasma particles. 

While it represents an important step forward, this description of the dynamics of 
heavy quarks by a Schrodinger equation and an effective potential has limitations. The 
potential is calculated, and well defined, only in the limit of infinitely massive quarks. 
Moreover, a simple potential description emerges only at large times, that is, at time scales 
that are large compared to the typical times characterizing the response of the plasma to 
perturbations. In the situations which we want eventually to deal with, namely the fate 
of bound states of heavy quarks in the environment created in ultra-relativistic heavy ion 
collisions, all relevant time scales are mixed (see for instance [17]), and a description of the 
dynamics beyond that provided by a simple Schrodinger equation is called for. This paper 
represents an attempt in this direction, building on the approach developed in [13]. Our 
strategy, already sketched in [18], is the following. The heavy quarks are treated as massive, 
non relativistic, particles coupled to a fluctuating gauge field. The dynamics of the heavy 
quark is then encoded exactly in a quantum mechanical path integral, while the average 
over the gauge field fiuctuations is entirely determined by the properties of the medium. 
If one restricts oneselves to approximations where this average is Gaussian, and hence can 
be performed analytically, the gauge fields can be eliminated completely, leaving a path 
integral for a non relativistic particle with a non local (in space and time) self-interaction 
term. This path integral is reminiscent of that introduced by Feynman in his treatment 
of the "polaron" [19]. An approximation that leads to a Gaussian average (at least in the 
Abelian case), is the hard thermal loop approximation (HTL) [20]. We shall make use 
of such an approximation, because of its simplicity, and also because it encompasses the 
dominant plasma effects that one wants to include: screening effects, collective modes, and 
collisions. We emphasize, however, that this approximation, which concerns primarily the 
medium in which the heavy quark propagates, can be improved without altering the basic 
structure of the problem. 

The present paper has an exploratory character and represents only a first step in this 
long-term goal. It focusses on the dynamics of a single heavy quark moving in a plasma of 
light charged particles, for which we provide a simple model. We use for the quark-gluon 
plasma an idealization where only Abelian (in fact. Coulomb) interactions are taken into 
account. We also assume, for simplicity, that the plasma particles are fermions, i.e., quarks. 
In short, we model the quark-gluon plasma by an electromagnetic plasma, treated within 
the HTL approximation. This is enough to take into account typical plasma effects, such 
as screening. Landau damping of collective excitations, and collisions between the heavy 
quark and the plasma particles. These phenomena are characterized by a single scale, 
the Debye screening mass m/j, to which, in our numerical studies, we shall give a value 
characteristic of a quark-gluon plasma at a given temperature (thereby taking effectively 
gluons into account). The dynamics of the heavy quarks is then treated exactly within a 



path integral of the type discussed above, with a non-local self-interaction whose space-time 
properties are controlled by the Debye mass. 

Our paper is organized as follows. In Sect. 2 we establish the general setting: the basic 
properties of the propagator of a heavy particle are recalled, a description of the medium of 
light charged particles in which the propagation takes place is given, the path integral for 
the heavy quark propagator is constructed. This path integral is calculated exactly in the 
limit of an infinitely massive quark, and known results are recovered concerning the long 
time behavior of the heavy quark propagator in this limit. Then, in Sect. 3, we calculate 
the heavy quark propagator in the one-loop approximation, providing a detailed analysis of 
the coupling of the heavy quark to the collective plasma excitations and of the role of the 
collisions. We also calculate the spectral function and the resulting Euclidean correlator. In 
Sect. 4 we present the results of the numerical evaluation of the path integral in Euclidean 
space-time, using Monte Carlo (MC) techniques. We use the Maximum Entropy Method 
(MEM) to reconstruct the spectral density. Within our present implementation of this 
method, we can only reproduce, semi-quantitatively, the main features of the spectral 
density. Finally, Sect. 5 summarizes the conclusions. In Appendix A we give a self- 
contained presentation of an exactly solvable toy-model which captures general features of 
the heavy quark correlator and its spectral function, and this for any value of the coupling 
constant.^ Appendix B provides details on the numerical evaluation of the path integral. 

2. A path integral for the heavy-quark propagator 

In this section, we recall general properties of the heavy quark propagator, and establish 
the basic path integral that describes the dynamics of the heavy quark coupled to a gauge 
field that is integrated out via a Gaussian averaging. 

2.1 The heavy quark-propagator. Generalities 

Most of the physical information that we are interested in can be obtained from the study 
of the following correlator 

G>{t,r\0,O) = m,r)iP\0,O)), (2.1) 

where 'ip{t,r) denotes the heavy quark field. In the following we shall most of the time 
use the simplified notation G^{t,r) for G^{t,r\0,O). The expectation value in the above 
formula is a thermal average over the states of light particles (with no heavy quark present) 
that will be specified later. At this stage, we simply note that the full Hamiltonian H can 
be decomposed into three contributions: 

H = Hq + Hmf.fl + Hint, (2.2) 

where Hq is the (non relativistic) Hamiltonian describing the heavy quark in vacuum, 
Hfned is the Hamiltonian of the medium in which the heavy quark propagates, and Hint 



^A somewhat similar model, with however different emphasis, was considered in Ref. [21]. 



represents the interactions between the medium and the heavy quarks. For the parts that 
depend exphcitly on the fermion field, we have 

Hq = M j d'r^\r)^{r) + j d'r^\r){-^\^{r), (2.3) 

and 

Hint =9 / d'ri;\r)i,{r)AQ{r), (2.4) 



where AQ{r) represents the local electrostatic potential created by the light particles. The 
full Hamiltonian commutes with the number of heavy quarks Nq: 

[H,Nq]=Q, Nq= f Sri,\r)i,{r), (2.5) 



and one can classify its eigenstates according to the number of heavy quarks that they 
contain. In particular, one may write a spectral decomposition of the correlator (2.1): 

G>(t,r) = V^^e*(^"-^'^)*(n|V'(r)|m)(m|V'^(0)|n), (2.6) 

n,m 

where the states \n) contain no heavy quark, while the states \m.) contain one heavy quark, 
i.e., 

^(r)|n) = 0, NqItti) = \m). (2.7) 

In eq. (2.6) Z is the partition function of the system without heavy quark. It follows also 
from Eq. (2.5) that G'^{t,r) = {il;^ {0,0)ip {t,r)) = 0, so that the time-ordered propagator, 
G{t) = i {T^{t)i>^0)) = ie{t)G>{t) -ie{-t)G<{t), and the retarded propagator Gi?(t) = 
ie{t)[G>{t) + G<{t)] are identical, G{t) = GR{t) = ie{t)G>{t). Similarly, the spectral 
density is given by the Fourier transform of Eq. (2.1), namely: 

/oo 
(ite^^*G>(t), (2.8) 

-oo 

with the inverse relation 

/ + 00 J 
^e-'^V(w). (2.9) 

-oo 



27r 
In this last equation we have exploited the analyticity of G^{t) in the strip —f3 < Imt < 0, 



and set t = —it, with < r < /3. Inverting this relation, namely calculating cr{uj) from 
the Euclidean correlator G^{—iT) is a difficult (well known) problem that we shall address 
briefly in the last part of this paper. 

By noting that Hmed does not depend on ^, one finds (with all fields in the Heisenberg 
representation) 

[^, H] = [^, Hq + Hint] = (^ - ^ + 9Ao{t, r)\ V^(t, r), (2.10) 

so that, from the equation of motion idtip{t,r) = [ip,!!], we get 

idtG>{t,r) = (m - ^^ G>{t,r)+g{Ao{t,r)i;{t,r)4^Ho)). (2.11) 



In the absence of interactions, this equation has the famihar solution 



M \ 



3/2 



G>(t,r) = e--'^^J e-^-'/^*, (2.12) 

corresponding to the initial condition Gq^I = 0, r) = 5{r). Note that this initial condition 
still holds in the presence of interactions, i.e., G^{t = 0, r) = 5{r), as is easily verified. A 
further exact relation is obtained by considering the equation (2.11) at t = 0: 

i dtG>it, r)l^^ =(^M-^ + g{Aoit = 0,r))) 6{r). (2.13) 

Since the thermal average involves only states of the medium which do not contain heavy 
quarks that could polarize it, we have (ylo(r)) = (1/Z)^„e~'^ "(n|ylo(r)|n) = 0: the 
interactions do not contribute to the leading (linear) order in a small time expansion. 
Pushing this expansion to second order, one gets 



-d^G>it,r)l^^ 
or, taking a Fourier transform 



2X2 



2m) 



Sir), (2.14) 



-a/G>(t,p)|^^,= ^M + ^J +g'{At,), (2.15) 

where (^q) stands for {AQ(t = 0,r = 0)). Thus at order t^, the effect of the interaction is 
governed by the size of the fluctuations of Aq, an intrinsic property of the medium to be 
discussed further later. Note also that the coefficient of t^ in the expansion of G^(t,p) at 
small t is of order g^ . 

We may turn these relations for the derivatives of G^{t,p) at t = into sum rules 
for the spectral density. From the initial condition (see after Eq. (2.12)), and the relations 
(2.13) and (2.14) above, one gets, respectively, 

doj , , f^ duj ^ ^ ^^ p"^ 

^^cx(a;,p) = l, y_^--(-^P) = M + ^> 

£^.M^,p)^(m + ^)%.^(^^ (2.16) 

The last sum rule assumes that {Aq) is well defined. However, as we shall see in the next 
subsection, within the approximation used in the present paper (Aq) is in fact given by a 
divergent integral, so that this sum rule will not apply. Accordingly the short time behavior 
of the correlator will not have a simple Taylor expansion as assumed in the discussion above 
(Eq. (2.15)). 

The Fourier transform of G^(t, r) used above (see Eq. (2.15)) is of the form 

G>{t,p) = Jd\e-'P'-G>{t,r) = (0^(^)0^), (2.17) 



where Op and Op are the Fourier transform of the field operators ^(r) and ip'{r), respec- 
tively, and we have used the translation invariance of the medium in order to implement 
the conservation of the total momentum ((apoL) ~ ^p,p')- For the value t = — i/3, the 
correlator G-^{—i/3,p) yields the difference of the free energies of the systems with and 
without a heavy quark. To see that, note that this free energy difference is given by 

exp[-/3AFQ,p] = —^{n\ape''^^al\n) 

n 

= ^5;e-'^^"(n|ap(/3)at,(0)|n) = G>(-i/3,p). (2.18) 

n 

In the first line of Eq. (2.18), the states ap\n), while not eigenstates of H, constitute a 
basis of states with momentum p and containing one heavy quark. Thus, the sum over the 
states \n) in the first line of Eq. (2.18) is indeed the partition function for the system with 
one heavy quark and total momentum p. 

In preparation for the forthcoming discussion, let us recall that the propagator of the 
heavy quark, treated as a non relativistic quantum mechanical particle, may be given a 
path integral representation [19]. With AQ{t,x) considered as a given external potential, 
we can write (t > 0): 



G>{t,r)= / Vz exp i dt' (-Mz^ - gAo{t,zU 



(2.19) 



where the symbol /q Dz indicates a path integration over paths z(t) such that z(0) =0 
and z{t) = r. The transcription of this expression in imaginary time reads (r > 0): 



G>(-ir,r)= / Vzexp - f dr' (^Mz'^ + igA^{T,z)j 



(2.20) 



where, aside from making the familiar substitution t — t- — ir, we have also introduced the 
Euclidean field Aq^TjV) = —iAo{t = —iT,r). 

2.2 A model for the medium 

The medium is modeled by a plasma of light fermions with Coulomb interactions. Because 
of its large mass, the heavy quark has a small velocity, and consequently its ability to induce 
magnetic excitations of the medium is small; accordingly these magnetic excitations will 
be ignored. The Hamiltonian reads then 

H^ed = fd\(\r)hoar) + I fd\d\'p{r) f p{r'), (2.21) 

J 2 J 47r|r — r'\ 

where ^(r) and ^^(r) denote the field operators of the light fermions, p{r) = £,'^{r)S,{r) is 
the charge density of the light particles, and /iq their free Hamiltonian. 

To the full Hamiltonian of the system corresponds an Euclidean action of the form 

S = Sq + Sint + Smed-, with 

Smed = I C {X) {dr + Hq) ({x) + ^ [ p{x)K{x, x')p{x'), (2.22) 



9 

X ^ Jx,x' 



and 

Srnt = 9^ I Pq{x)K{x,x')p{x'). (2.23) 



x,x' 



We have set Pq{x) = ip*{T,r)7p{T,r), and 

d*x= I dr d^r, x={T,r). (2.24) 

The operator K{x,x') is given by 

-VlKix, x') = 6{x - x'), K{x, x') = 5{t - t') —-^ -. (2.25) 

47r|r — r'\ 

In calculating the partition function of the system, one can proceed in a familiar way, and 
integrate over the light fermions after eliminating their density p{x) in favor of a gauge 
potential Aq{x) {—igK ■ {p + pq) — )■ A^). One then obtains 

fv{C,0(^~^^'"'^''^"""'^= f VA^ e^''^^^o] ^ (2.26) 

where 

S[Ai] = ig [a^{x)pq{x) 



X 



- Trln{dr + ho + igA^) + 1 f Ai{x)K-\x,x')A^{x'), (2.27) 

^ Jx,x' 

and the field A^ obeys periodic boundary conditions in imaginary time, yl^(0,r) = 
Aq^/S, r), reflecting the fact that the medium of light particles is in thermal equilibrium at 
temperature T = 1//3. 

At this point we perform the main approximations that will yield a simple model for 
the medium. In the expansion of the fermionic determinant (the second term in the r.h.s. of 
Eq. (2.27)) in powers of A^ , we keep only the quadratic term. Furthermore, we keep only 
the leading high temperature approximation to the corresponding 2-point function (the 
so-called hard-thermal-loop (HTL) approximation [20]). Note that in the case of QED, the 
HTL approximation automatically truncates the expansion of the determinant at quadratic 
order. Further discussion of the validity of this approximation will be made shortly. At 
this point we note that once it is done, we are left with a simple quadratic action: 

S[Ai] =tg [ A^{x)pq{x) + \j A^{x)K-\x,x')A^{x'). (2.28) 

J X ■^ J x,x' 

The propagator A(x,x') = {Aq{x)Aq{x')) is given in Fourier space by K^'^{z,q) = q^ + 
11(2:, g), where Ii{z,q) is the (longitudinal) polarization tensor in the Coulomb gauge: 

n(z, q) = ml{l- Q{z/q)) , Q{x) ^ ^ In ^, (2.29) 

2 X — I 

with rriD = II(z = 0, q) is the Debye mass. The Debye mass is the mass scale that 
characterizes the response of the medium. 
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Figure 1: The function A(r, r) as a function of r/rn — rmu for different values of t//3 (decreasing 
from bottom to top ). Note that as long as t ^ 0, A(t, r == 0) is finite. However, A(0, r) diverges 
logarithmically as r — >■ 0. 



At this point, we note that the equations we have written hold exactly only for the 
hot electromagnetic plasma. However, at this level of approximation, the main difference 
with a quark-gluon plasma lies in the value of the Debye mass that, in a QCD plasma, 
receives also contributions from gluons. In the numerical studies to be presented below, in 
order to get orders of magnitudes that are relevant for the quark-gluon plasma, we shall 
adjust the Debye mass to the value it would have in a quark-gluon plasma at the considered 
temperature, that is we shall use the QCD HTL expression tud = qIT'^ (-^c/3 + Nf/6), 
with gl/^'K = as the strong coupling constant. With ag = 0.3, Nc = 2> and Nf = 2, this 
yields a value mo = 713 MeV for T = 300 MeV. Furthermore, the coupling of the heavy 
quark to the plasma particles involves qI/^t^ multiplied by the Casimir factor Cp = 4/3. 
We shall absorb this factor Cp into the coupling g, denoting the product glCp/'^T^ by 
a = g"^ /A.TT. Thus a coupling constant a = 0.4 in our notation, corresponds effectively to 
Us = 0.3 in QCD. 

The propagator A~^(z,q) introduced above contains all the information about the 
screening phenomena and the response of the medium to the presence of the heavy quark. 
It differs by a sign from the longitudinal gluon propagator in the HTL approximation 
(called ^L in Ref. [13]). It is convenient to subtract from the latter the instantaneous 
Coulomb interaction which would contribute here only to the self-interaction of the heavy 
quark. Thus we define 

1 1 

q2+n(z,q) ~^ 



A(z,q) 



(2.30) 



This new object A(z,q) is proportional to xiz,q), the density-density correlation function 
of the medium: A{z,q) = {l/q'^)x{z,q). One has, in a mixed representation: 



A(T,g) = 1^" |^e-5o-p^(go,g)[0(r) + Niq% (2.31) 

where the spectral function pL{u},q) reads [13] 

PLioJ, q) = 2tt {ZL{q) [5{u-u:L{q)) - 5{uj+u:L{q))] + 9 {q^ - uj'') p l{uj , q)] . (2.32) 

It displays two types of contributions: A continuum term arising from the imaginary part 
developed by the logarithm in Eq. (2.29) for space-like momenta, and which corresponds 
physically to scattering processes, and a pole term, coming from the solution, for time-like 
momenta, of 

q2+n(c^L(g),q) = 0, (2.33) 

which corresponds to an undamped plasma oscillation. Note that the residue Z^lq) quickly 
dies out as q grows beyond mi): 

ZL{q) ~ 3^ exp - , ^ ) , g > mi5. (2.34) 



Collective modes exist only for q ^ rao. 

The approximation that we are using to describe the medium to which the heavy 
quark is coupled is motivated by its simplicity, and also by the fact that it encompasses the 
important physical phenomena that characterizes weakly coupled plasmas, and that one 
wants to take into account: polarization and screening effects, collisions with the plasma 
particles. The latter, however, are not treated properly in the HTL approximation, and 
this will introduce (small) unphysical features in our results. As a concrete illustration of 
the difficulty we are referring to, consider the function A(r, r) that will play a central role 
in our calculations. This function can be obtained by a Fourier transform of Eq. (2.31) over 
spatial momentum, and it is displayed in Fig. 1. As indicated in the caption of this figure, 
A(0,r) is logarithmically divergent as r — )■ 0. This divergence is that of the fluctuation 
(^q) = A(0, 0), given by the integral 

{AD = j ^^PL{q\q)N{q^), (2.35) 

and comes form the continuum part of the spectral function (the contribution of the plas- 
mon is finite, due to the vanishing of the residue for large wave-numbers; see Eq. (2.34)). 
As already mentioned, the continuum part of the spectral function describes scattering 
processes, involving space-like gluons with energy w, momentum q. In the HTL approxi- 
mation, the phase space for these processes is given by |a;| < q (see Eq. (2.32)), i.e., it grows 
without bound as q increases, leading eventually to a divergence. An analogous divergence 
also occurs in the pair correlation function at short distance when this is calculated using 
the Vlasov equation (which is equivalent to the HTL approximation [22]). This is a well 
known difficulty in plasma physics (see e.g. [23]), and it can be cured by a better treatment 
of the collisions involving large momentum transfer. Indeed, the HTL approximation is 
valid only in the regime where 9 <C p, where p ~ T is a typical loop momentum (i.e., the 
typical momentum of the colliding plasma particles). A proper treatment of the collisions 



with g ~ p would lead to a modified phase space and a finite value for (Aq). For instance, 
in a full one-loop calculation, the phase space is given hy —q < u < q for q ^ p, but 
q — 2p < oj < q for q > p. A possible way to improve the calculation would be to introduce 
a cut-off to separate soft and hard contributions, and apply in each sector appropriate ap- 
proximations. We shall not do so here, because A{t,r) enters the calculation of the heavy 
quark correlator only through an integral so that its logarithmic singularity is tamed, and 
its physical consequences mild. We note however that the divergence of {Aq) modifies the 
small r behavior of the heavy quark propagator, and in particular it invalidates the Taylor 
expansion dicussed at the end of Sect. 2.1, beyond the linear order. 

2.3 Path integral representation 

We are now in position to write the propagator of the heavy quark in the form of a path 
integral. Gathering the results of the first two sections, we can write 

G>i-iT,r) = fvA^ exp \-\ f Ai{x)A-'ix,x')A^ix'] 

J L ^ Jx,x' 



Vx exp 





/ dr' [\mx^ + ieA^ (r, x) J 



(2.36) 



This path integral summarizes the model that we use. The dynamics of the heavy quark 
in a hot plasma is that of a non relativistic quantum particle moving in a fluctuating 
field Aq, and this is treated exactly by the path integral. The approximations only enter 
the description of the fluctuations of the held Aq which we assume to be Gaussian and, 
as we have just discussed, dominated by long wavelengths and low frequencies. Thus, 
any improvement of the description of the medium will affect only the first part of the 
functional integral (2.36), that is the weight of the integration over the field ^O; but it will 
leave intact the second part describing the motion of the heavy quark in the fiuctuating 
field. This is an important feature of the present description. 

As we mentioned earlier, it is convenient to subtract from the correlator A the con- 
tribution of the Coulomb interaction. This is most easily done after having performed the 

E 
' 



Gaussian integral over Aq , whence we can just replace A by —A. One gets 



G>{-iT, r) = r Vz e-^[^'^ , (2.37) 

where S'[2;,r] = 5*0 [2;, r] — F[z,t], with 

So[z,t]= r dT'\Mz\ (2.38) 



'0 

and 

^2 



F\z, t\ = %- r dr' r dT"A{T' - t", z(t') - z{t")). (2.39) 

2 Jo h 

The real time version of this path integral is obtained by replacing r by it, and substituting 

-5[z, r] in Eq. 2.37 by iS\z, t\ = i{So + F) with 

F[z, t] = ^ [ dt' [ dt"D{t' - t", z{t') - z{t")), (2.40) 

2 Jo Jo 
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and we have used [13] 

A(r = it, r) = -iD{t, r). (2.41) 

The correlator G^{t, r), when expressed in terms of the dimensionless variables tT and 
rT is of the form G^(t,r) = T^f{M/T,mD/T,tT,rT), with / a dimensionless function. 
When niu — )■ this reduces to the free propagator. Note that at fixed value of the coupling, 
771/5 /T is fixed, and G^ (t, r), when t and r are expressed in units of the inverse temperature, 
depends only on the ratio M/T. We shall refer to this scaling property of the correlator 
repeatedly. 

The parameter M/T controls the "diffusion", described by the correlator (2.12) in 
imaginary time: the smaller M/T, the more the heavy quark will move away form the 
origin in a given time. Note that this diffusion inhibits the effects of the interaction: 
because A(r, z) < A(r, 0) (see Fig. 1), the interaction favors paths for which z remains 
small (their weight in Eq. (2.37) is largest). 

One may also understand the effect of the interaction in the following way. The heavy 
quark produces a polarization cloud of light particles around itself. This induced charge 
screens that of the heavy quark over a distance scale of order m^ . When the heavy quark 
moves, its polarization cloud tries to adjust and follow its motion, but this takes time 
(of order m^ ). The heavy quarks sees then a restoring force produced by the lagging 
polarization cloud, which limits its excursion. 

In the limit M/T — t- oo, studied in detail in the next subsection, the heavy quark 
is frozen at it initial location: there is then no diffusion, and the effect of interactions is 
maximal. 

2.4 The infinite mass limit 

When M — 7- oo, the path integral can be calculated exactly. This is because, in this limit, 
the motion of the heavy quark is frozen and F becomes independent of the coordinates. 
Thus, in the infinite mass limit, the heavy quark correlator takes the form 

G> {t, r) = 6{r) e"^^* e*^W (2.42) 

where the function F{t) is the functional (2.40) restricted to z = 0: 

F{t) = %- f dt' f dt"D{t' - t",0). (2.43) 

2 Jo Jo 

The factor exp {iF{t)) mulitplying in Eq. (2.42) the infinite mass limit of the free correlator 
(2.12), summarizes the effect of the interactions. One can express F{t) in terms of the 
Fourier transform D{uj,q) of the (time-ordered) gluon propagator [13]: 

dq° pL{q°,q 



One gets 



»'"■■'> -/^ /-t+l) ^'^^'"'"'^'"*^ <'^''' 



^/ N n I duj l — cosiut) f d^q ^. , , ^^, 

m =9^1^ -^^ J J^ Di., q). (2.45) 
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It follows that at short times 

^2 



nt) ^\e j^j j^Dioj, q) = '-t'D{t = 0, r = 0). (2.46) 



For large time we use 



to obtain 



2 
J 


f duj f d?q 9'^^2r,f. 




l-cos(wt) 

lim 7: = TTldlUj). 


n 


2 

t) ^ ^tD{oj = 0,r = 0) = -tVop 



(2.47) 



(2.48) 

An alternative way to obtain this result is to start directly from Eq. (2.43) and to change 
variables t' — t" — )■ u, {t' + t")/2 — )■ T, and to observe that at large time t, one may integrate 
freely over u, thereby filtering out the zero frequency part of D{uj,q). This yields again 
Eq. (2.48). Thus, the large time {t ^ m~jj ) behavior of the system is determined by the 
static (oj = 0) response of the medium. Since, at large times, F is linear in time, Eq. (2.11) 
for G^{t,r = 0) is a closed equation that takes the form of a Schrodinger equation [13], 
with an "optical potential" Vopt given by 



9 f dq 



1 1 . Trm^T 



q'^+mj-, q^ IqKq^ + "I-d)^ 



2 2 

where we have used the susceptibility sum rule 



2 J (27r)3 
^m^ _ i^^ (2.49) 



f°° dqo PLiqo,q) ^ mjj 

J -00 27r qo q^{q^+ml) 



D> 

in order to perform the go integral needed to calculate Diyj = 0,g). Thus, the optical 
potential involves a real correction to the mass of the heavy quark, and an imaginary part 
that takes into account the coupling of the heavy quark to the complex configurations of 
the medium. Alternatively, one may view this imaginary part as reflecting the collisions 
of the heavy quark with the particles of the medium. As we shall see in the next section, 
T4pt can be identified with the one-loop on-shell self energy in the infinite mass limit. 

This imaginary part does not appear in the Euclidean correlator calculated at r = /3, 
which exhibits only the mass shift [13]: 

-rinG>(t = -i/3,p) = M--raj:) = M - — A(iw„ = 0,r = 0). (2.51) 

Note that in the present, infinite mass, limit, G^(t = — i/3,p) in fact does not depend on p. 
More generally, the Euclidean correlator has the form (2.42) with i-F(t) replaced by F(t), 
with (see Eq. (2.39)): 

F(r) = ^ r dr' f dr" ^{t' - r", 0). (2.52) 
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Since the dominant effect of the interactions can be characterized by the free energy shift 
(2.51), it is convenient to separate the corresponding linear growth of F{t), and write 



with 



and 



F(r)=Fi(r) + F2(r), (2.53) 

v( \ - s"^ f J^lJL f "^^o ^-^(^0' ^) cosh(go(T - (3/2)) - cosh(/3go/2) 

^'^^> - T y (2^ y ^^e smH(3qo/2) " ^^■^^> 

The function F2{t) vanishes at r = and r = /3. It is symmetric around r = /3/2, a 
property that follows immediately from the fact that A(r, q) depends only on |r|, and is 
periodic, A(/3, q) = A(0, g). Note also that the slope of F2{t) at the origin is equal and 
opposite to that of -Fi(''"), since that of F{t) vanishes. 

Many of (but not all) the features of the present M — t- oo limit are shared by the toy 
model presented in Appendix A, where one can find a more extended discussion of some 
of the points addressed in this subsection. 

3. One-loop calculation 

In this section, we present the results of the one-loop calculation of the heavy quark cor- 
relator. This provides insight into the dynamics of the heavy quark when the interaction 
is weak enough for perturbation theory to be applicable. All the numerical results to be 
presented are obtained with the value a = 0.4 of the coupling constant, which appears to be 
a moderate value for which the one-loop approximation remains reasonably accurate. This 
calculation, together with the exact large M limit that we have just discussed, will serve as 
a reference when discussing the results of the Monte Carlo evaluation of the heavy-quark 
correlator in the next section. 

The one-loop calculation is easier in momentum space than in coordinate space. To 
proceed we consider the analytic propagator 

where Ep = M + p'^ /2M, and the one-loop self-energy S(z,p) is given by the diagram 
displayed in Fig. 2. The retarded propagator is obtained as usual by setting z = uj+irj, with 
CO real. The imaginary part of the retarded propagator yields the heavy-quark spectral 
function 

a(uj,p) = 2lmG^(uj,p) = -, ^ ^. i^ — 7^777 ttt, (3.2) 

^ ^^ ^ '^' [tJ-Ep-ReS(w,p)]2 + r2(a;,p)/4' ^ ' 

where r(a;,p) = — 2ImS'^(a;,p) = — 2ImS(z = a;-)-i?7,p). Eventually the Euclidean correlator 
will be calculated using Eq. (2.9). Since we shall consider only the case p = 0, we shall use 
in most of this section the simplified notation S(z) for S(z,p = 0), and similarly for other 
related functions. The relations (3.1) and (3.2) are general, but in the rest of this section, 
G and S will refer to one-loop quantities (unless stated otherwise). 
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Figure 2: The one- loop self-energy diagram for the heavy quark. The blob on the interaction line 
reminds that the latter represents a resummed HTL propagator of a longitudinal gluon. 



3.1 The one- loop self-energy 

A standard calculation, implementing approximations that are valid when T/M <^ 1, yields 
the analytic one-loop self-energy 



S(^) = g' 



dk 



-pL{k ,k)- 



2tt 



Ek-kO- 



(3.3) 



Expressing momenta and energies in units of T, one sees that T,{z) is a function of the 
form T,(z) = Tf{z/T,M/T,mD/T). At fixed value of the coupling constant, mo/T is 
a constant, so that, the only relevant parameter is the ratio T/M, as we have already 
mentioned. 

By using the explicit expression for the gluon spectral function pL{k^,k) given in 
Eq. (2.32), one can re- write Eq. (3.3) as 



S(^) = / 



dk 



zUk) 



l + N{iUL{k)) 



+ 



NiuLik)) 



+ 1 ^27r h{k',k) 



ujiik) z-Ek+ujL{k) 
1+N{k^) 



+ 



+ 



N{k^) 



(3.4) 



z- Ek 
'^ dk^ 

Ek-k^ z-Ek + k^_ 

This expression exhibits two types of contributions that are are illustrated in Fig. 3: a pole 
contribution whose energy denominators are associated with processes of emission or ab- 
sorption of collective plasmons by the heavy quark, and a continuum contribution coming 
from the continuum part of the gluon spectral density; the latter contribution represents 
the effect of collisions between the heavy quark and the particles of the medium, medi- 
ated by space-like gluons. It is convenient to evaluate separately these two contributions. 
Accordingly we set r(a;) = rP°'^(a;) + r'=°"t(u;). 
For the pole contribution one gets: 



■ppole/ 



OJ 



9 



TT 



kl 



\K+^'L{kl)\ 



ZL{ki)[l + N{uL{ki))] + 



+ 



^ \FJ 



^2 



UJ 



Xk2)\ 



ZLik2)N{0JLik2)) 



y (3.5) 
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Figure 3: The different processes contributing to the imaginary-part of the heavy-quark seff- 
energy: (a)-(b) emission-absorption of a plasnion and (c)-(d) collisions with the plasma particles, 
mediated by one-gluon exchange. 



where fci and /c2 are implicit functions of uj given by 



U! 



Eki +a;L(/ci), 



OJ 



Ek2 - Wl(/C2 



(3.6) 



and the primes in the denominators of Eq. (3.5) denote derivatives with respect to k. 
Here ujL{k) is the plasmon dispersion relation (see Eq. (2.33)), whose behavior for small 
momenta reads: 



UJL{k) ~ w 1 



^l'' 



3 k'^ 

ujiik) ~ Wpi-F— — . 



(3.7) 



The solutions of Eqs. (3.6) can be read out from Fig. 4 where the two curves Ek^LOL{k) are 
plotted as a function of k. The first equation (3.6), oj = E^ + 0JL{k), has a single solution 
starting from the plasmon-emission threshold uj = M + ujp\. The number of solutions of 
the second equation depends on the ratio ojpi/M. From Eqs. (3.6) and (3.7) one sees that 
for Wpi > |M the dispersion relation lo = Ek — ujL{k) starts with positive curvature and 
it contributes to FP"^*^ with a single solution starting from oo > M — Wpi. In the case of 
interest, T/M ^ 1, we have Wpi < |M, and there are two solutions for M/2 ^co < M — ujp\ 
and only one for io > M — ujpi. Note however that any contribution corresponding to 
fc ~ M is damped by the plasmon-residue. As a result, there is effectively no pole term 
contribution for M — ujp\ < oj < M + cOpi, as it can be seen from the left panel of Fig. 5. 
From Fig. 4 one realizes also that there are values of co for which the dispersion relations 
uj = EkiujL{k) display stationary points which, because of the denominators in Eq. (3.5), 
lead to singularities in F^ (oj) [Van- Hove singularities), clearly visible in the left panel of 
Fig. 5. Notice that, as the ratio T/M gets larger, the one occurring at u;r^M/2 acquires 
more importance, getting less suppression from the plasmon residue. 

The continuum contribution involves the spectral density /3l(/cO) k) which has support 
for \ko\ < k. It follows then from Eq. (3.4) that the continuum contribution to the imaginary 
part of S comes from values of k,uj such that —k < uj — M — ^^ < k. The boundaries of 
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Landau-damping 
plasmon T=300 MeV 
plasmon T=600 MeV 



k/M 



Figure 4: The dashed curves represent, for two different temperatures, the functions Ek±ujL(k), 

with Ek = M+^, and ^^(fc) the plasma dispersion relation. For fc = 0, a; = M±ujpi, with 

cjpi = 7Ti£)/v3 the plasma frequency, proportional to the temperature. Here, M — 1.5 GeV and 

Wpi = 412 MeV for T = 300 MeV. The two full lines delineate the support of the continuum part 

fc2 



of the gluon spectral function, that is the region —k < oj — M— ^ < k. The largest temperature 
r = 600 MeV corresponds to a plasma frequency a;pi = 824 MeV, very close to the "critical value" 
3/5M discussed in the text. 
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Figure 5: The pole (left) and continuum (right) contribution to T{lu) for a quark mass AI — 1.5 
GeV, and various temperatures. The plasma frequency is proportional to the temperature and has 
values Wpi = 275, 412 and 687 MeV for the temperatures T = 200, 300 and 500 MeV, respectively. 
Notice, in the pole contribution, the Van- Hove singularities and the gap, for M — uipi<oj<M + uipi, 
that increases with temperature. The continuum contribution grows linearly with temperature, and 
the threshold at cj = M/2 is clearly visible. 



this domain are displayed in Fig. 4, and r™'^*(a;) is given by 



ncont 



H 



9 



M+M{uj) 



e{uj-M) _ edkl3L{uj-Ek,k)[l+N{co-Ek)] + 

71" [ J~M+M{u) 

+ e{uj-M/2)e{M-uj) I _ k'^dkf3L{io-Ek,k)[l+N{uj-Ek)]}, (3. 

Jm-M{lo) \ 
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Figure 6: Imaginary part (left) and real part (right) of the self-energy S. The horizontal lines, 
labelled "static limit", indicate the values of r(M— )-oo) and ReE(-M— )-oo). With the parameters 
a = 0.4 and T = 300 Mev, these are respectively 120 MeV and -143 MeV. Within the gap ±ujpi, T is 
an increasing function of M, while ReS is a decreasing function of M. Both functions nearly reach 
the infinite mass limit when M = 45 GeV. 



where M{uj) = ^/AP+2M{lo-M). Note in particular the lower threshold at w = M/2, 
corresponding to the minimum at k = M of the lower boundary of the support displayed 
in Fig. 4. This is clearly visible in the plot of T^"^^{uj) in the right hand panel of Fig. 5. 

A similar analysis can be done for the real part of the self-energy. This will not be 
detailed here. We just present in Fig. 6 the result of the full calculation of the imaginary 
part (left panel) and the real part (right panel) of S, for different values of the heavy quark 
mass, including the limiting case of infinite mass. Let us recall that in the latter case, we 
have analytic results [13] for the on-shell values (corresponding to w = M). From Eq. (3.4) 
one sees that only the continuum part contributes to the on-shell imaginary part 



ImS^(M^oo) = -— lim / — -^iV(A;°)pL(fc°,fc) = ^ r(M^oo) = ar, 

2 ko^oj {2'ir)'^ 2 



while the real part receives contribution from both parts of ^(z): 



ReS(M^oo) = - — 



dk 

(2^ 



dk^ PLik^,k) 



amo 



27r 



fco 



(3.9) 



(3.10) 



These values, which coincide with the values obtained for the exact "optical potential" 
in Eq. (2.49), are indicated by the horizontal lines (labelled as "static result") in Fig. 6, 
while the curves representing the full expressions of ReS((x') and r(w) in the infinite mass 
limit are labelled as "M=infinite" . One sees from this figure that the infinite mass limit is 
nearly attained for M = 45 GeV, and that finite mass effects do not change the qualitative 
behavior of the self-energy. To get a quantitative measure of these finite mass effects, we 
determine the shift 5M = M' — M of the heavy quark mass as given by the solution of the 
equation 

M' - M = ReS(M'). (3.11) 
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Figure 7: One-loop spectral function as a function of w— M, for various values of the heavy quark 
mass, and a fixed tempearture T — 300 MeV. The curve corresponding to M = Ab GeV is hardly 
distinguishable from that representing the one-loop infinite mass limit. The smaller the mass M, 
the smaller the shift of the main peak. 



This can be obtained graphically, as the intersection of the line uj — M with S(cj) in Fig. 6. 
Values of the mass shift 5M thus obtained are reported in Table 3.1. 



T/M 





0.0067 


0.067 


0.133 


0.200 


0.333 


6M/T 


-0.407 


-0.4 


-0.357 


-0.335 


-0317 


-0.288 


^Fq/T 


-0.416 


-0.409 


-0.362 


-0.336 


-0.318 


-0.274 


ReS(M)/r 


-0.476 


-0.457 


-0.41 


-0.38 


-0.357 


-0.326 



Table 1: The mass shift 5M obtained from the solution of Eq. (3.11), the one- loop free energy 
shift AFq, and the real part of the on-shell self-energy ReI](Af) (which equals the exact energy 
shift in the infinite mass limit), as a function of T/M. 



One sees from this table that the larger the ratio T/M, the smaller the mass-shift. 
This is in line with what one expects from the effects of diffusion that increase as T/M 
increases, and inhibits the effect of the interaction. Note also that the mass shift obtained 
as the solution of Eq. (3.11) is numerically very close to the free energy shift calculated 
from the Euclidean correlator G{—i(3). It is smaller (in absolute value) than ReS(a; = M), 
as can be also directly seen in Fig. 6. 

3.2 One-loop spectral function and Euclidean correlator 

The spectral density can be readily calculated from the real and imaginary parts of the 
self-energy (see Eq. (3.2)). It is displayed in Fig. 7. The dominant feature is the existence of 
a main peak, approximately located at the value of a; = M', with M' given by Eq. (3.11), 
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Figure 8: The quantity -F'^^'^(t) (see Eq. 3.12)) for various values of the mass and fixed T = 300 
MeV. For the M — co case we also plot the function F{t) (the curve labelled "exp"), so as to get a 
measure of the accuracy of the one-loop approximation. The value of the free energy shift /S.F can 
be read on the left panel as the value of F'-^'^fS) and is reported in Table 3.1. In the right panel we 
plot the same quantities after subtracting the linear r-dependence driven by the free-energy, that 
is the function F^^'\t) in Eq. (3.13). 



as can be expected on general grounds from Eq. (3.2). In addition to the main peak, 
two secondary bumps appear in the spectrum at values of the energy a; ~ M =b Wpi, and 
come from the energy dependence of the imaginary part of S discussed in the previous 
subsection. The spectral density satisfies the sum rules (2.16): it is normalized to 1, and 
its first moment remains equal to M. Note that the infinite mass limit gives an accurate 
picture, only mildly modified by finite mass corrections, down to values of the mass of the 
order of 1.5 GeV. In particular finite T/M effects seem to be important mainly for the 
shift and the broadening of the main peak, affecting on the other hand very mildly the 
secondary bumps. 

By using the relation (2.9), one obtains from the one-loop spectral function the corre- 
sponding Euclidean correlator. This is displayed in Fig. 8, for different values of the heavy 
quark mass. What is plotted in Fig. 8 is actually the function 



F(iO(, 



In 



G>{-iT) 

Goi-ir) ■ 



(3.12) 



As it can be seen, all the curves, start with zero slope at r = 0. This is related to the 
general feature that the interactions do not introduce any corrections linear in r at small 
r, which in turn may be linked to the first two sum rules (2.16)) which are satisfied in 
the one-loop approximation. This represents actually an important consistency check of 
the numerical calculation, given the indirect way by which the Euclidean correlator was 
obtained. The value of the Euclidean correlator at r = /3 measures the free-energy shift AF 
caused by the addition of the heavy quark, and can be read off Fig. 8. As already obtained 
in the case of 5M, one finds a smaller shift as the ratio T/M gets larger (see Table 3.1). 

A different way to plot the Euclidean correlator is offered in the right panel of Fig. 8. 
There we have separated the linear r-dependence driven by the free-energy shift, writing 
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Figure 9: A comparison between the first-order perturbative expansion of the Euchdean correlator 
GpcrtlT) = G'o(t) + G'i(t) (dashed lines) and the full one- loop correlator G(t) obtained from the 
resummation of the Dyson series (continuous lines) . 



(see Eq. (2.53)) 



^(10(,) = F^^^\r)^F^^\r\ F^ir) = -^F^'Hf^)- 



(3.13) 



The difference of behavior that is observed is quite similar to that obtained in the toy 
model presented in Appendix A. Note in particular that the symmetry around /3/2 that is 
present in the exact M = oo limit is lost in the one-loop approximation (also in the infinite 
mass limit of the latter) . 

Finally it is of interest to study the accuracy of the vi^eak-coupling expansion at short 
time. To that aim, we expand the propagator to order a 

G(r) = Go(T) + Gi(r) + ..., (3.14) 

with Go(t) = e"*^"^, and Gi(r) is given by the one-loop self-energy: 

Gi{t) = g\-^'^ Tdr' r dr" / ^A(r' - r", fc) e-('=V2M)(r'-r"), (315) 

Jo Jo J (271") 

In order to calculate the time intergral, one may express the gluon propagator in terms of 
its spectral density. One gets then: 



Gi{r)/e 



^Mr ^ g2 



dk 



dk^ PL{k^,k)[l + N{k^)] 



(27r)3y 2tt k^ + k^/2M 



dk 



dk'^ PLik^,k)[l + N{k'^)] 
(27r)3 J 27r (A;0 + k^ /2Mf 



1 _ g-(A:"+fe2/2M)T 



(3.16) 
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The result, for the zero-momentum case, is plotted in Fig. 9 and compared with the full 
one-loop calculation. As it can be see, for the moderate coupling a = 0.4 considered here, 
the weak coupling expansion is accurate till quite large values of t//3. What is perhaps 
surprising is the dependence on the mass M, which reflects a non analytic behavior at 
small r. Assume indeed that a Taylor expansion of Eq. (3.15) exists. Then, the leading 
term in this expansion, of order r^, is obtained by setting r' = r" = in the integrand, 
leading to the result 

2 2 /. jr, 2 2 

Gi(r)/e-^- =^ / -7^A(0,fc) = ^A(r = 0,r = 0), (3.17) 



(27r)2 

which would be independent of the mass M. However, as already stressed, A(t = 0, r = 0) 
is divergent, so that Eq. (3.15) has no Taylor expansion. The integral over r' and r" in 
Eq. (3.15) exists however, and because of the exponential factor, it acquires a dependence 
on the mass M: it is largest in the limit M — )■ oo, and decreases as M/T decreases. This 
is the trend seen in Fig. (9). 

4. Numerical results: MC simulations and MEM analysis 

In this section we present the results of the numerical evaluation of the path integral for 
the heavy-quark correlator. We shall also discuss the spectral density obtained from the 
latter through an analysis based on the Maximum Entropy Method (MEM) [24]. Since no 
ambiguity can arise, we use in this section the simplified notation G{t, r) for the Euclidean 
correlator in place of G'^{—iT, r) used in the rest of the paper. This correlator is obtained 
from the path integral derived in Sect. 2.3. By taking the ratio of G{t, r) with the free 
propagator Go(r, r) (see Eq. (2.12)) one obtains 

Gir,r) _Io^ze~^ol^U^(^) _^^^,^^^ 



(e^^-'-J), (4.1) 



with 



and 



Go(r,r) /;pze-^oW 

So[z,t] = J^t'^Mz^ (4.2) 



F[z, r] = ^ r dr' r dT"A{T' - r" , z{t') - z{t")). (4.3) 

^ Jo Jo 

The functional F[z, r] is a known functional of the path, with A(t, r) an intrinsic property 
of the plasma, calculated as indicated in Sect. 2.2. The calculation of G{T,r) according to 
Eq. (4.1) amounts to an average that can be performed using Monte Carlo (MC) techniques. 

4.1 Monte Carlo evaluation of the path integral 

In fact, to proceed with the MC calculation, we shall take a slightly different route than that 
suggested by Eq. (4.1). This is because we want to include the effects of the interaction in 
the samples of paths used in the averaging. While this may not be necessary in the present 
one particle problem, this is essential when dealing with the two particle problem that we 
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Figure 10: Left panel: The quantity F^'^ij^r = 0) for various temperatures (and Af = 7.5). As 
the ratio TjM decreases the curves move closer to the static result. Right panel: F^'^ij^r') for 
T= 1 and various values of r. 



plan to address in the future. Thus, using a standard strategy, we define a propagator 
Ga{T,r) as in Eq. (4.1) but with the action replaced by 



Sa[z,T] = So[z,t] - aF[z,T], 



(4.4) 



with F[z,t] given by Eq. (4.3). Clearly, 5a[z,r] interpolates between 5o[2;,r], correspond- 
ing to a = 0, and the full action S'o[2;,r] — F[2;,r] reached for a = 1. By taking the 
derivative of In Ga with respect to a one obtains 



1 dGa (r, r) _JVz F[z] exp [-5, [z]] 



Ga{T,r) da fVzexjp[—Sa[z\] 

and G{t, r) is recovered after integration over a: 

In-^^ =lnf(e^[-M = ['da ^'^^ (r,r) 
Go{T,r) V ^ Jo da 



{nz])a, 



(4.5) 



da {F[z])^. 



(4.6) 



The a-dependent average value appearing in the right-hand side of Eq. (4.5) is evaluated 
using a MC algorithm whose details are given in Appendix B. 

4.2 The Euclidean correlator 

The heavy quark correlator is calculated first in coordinate space, and then at zero spatial 
momentum. Calculations have been performed for a fixed mass M = 7.5 and temperatures 
ranging from T = 0.75 to T = 2. Recall that all energies in the MC calculation are expressed 
in units of 197.3 MeV (so that T = 1 corresponds to T ~ 200 MeV, and M = 7.5 to 
M ~ 1.5 GeV). It is also useful to remember in the following that the ratio of propagators 
in the left hand side of Eq. (4.6) is a dimensionless function of t//3, r//3, T/M, and mu/T. 
Actually, since we keep the coupling constant a = 0.4 fixed, T/M is the only relevant control 
parameter. 
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For the ease of presentation we define 



pMC/ „\ _ ,„ G{T,r) 



F--(r,r) = ln-^^^. (4.7) 

Go (r,r) 

Tlie quantity F^*-'(t, r) is displayed in Fig. 10 together with its infinite mass limit, the 
function F{t) (see Eq. (2.52)). Since A(r, z = 0) > A(r, z), F[z,t] < F{t): hence diffusion 
tends to decrease the magnitude of F^^(t, r). Thus, the larger the ratio T/M, the larger 
the diffusion, and the lower is the corresponding curve in the left panel of Fig. 10. The 
panel on the right hand side of Fig. 10 indicates that the effect of the interactions depends 
mildly on r: it attenuates very slowly as r increases. 

We now consider the correlator projected to zero momentum 

GiT,p = 0)= f drG{T,r). (4.8) 

It is again convenient to study the ratio 



^j^MC(,,p=o) ^ G{t,p = 0) ^ J dr exp[-MrV2T] Gjr, r)/GoiT, r) 
° ^Go(r,p = 0) /drexp[-Mr2/2T] 



(4.9) 



This expression lends itself to a convenient numerical evaluation. Indeed, as an outcome 
of the MC simulations, for each r, one knows G/Go(t, r) for a discrete, and rotationally 
symmetric, set of values {vi} (typically rj^2 fm). One can then write 



G(r, p = 0) Z^ rf exp[-Mr2/2r] G(t, r,)/Go(r, r, 



Go(r, p = 0) Zi rj exp[-Mr2/2r] 



(4.10) 



which is the formula used to obtain G(r, p = 0). 

A further remark is in order. As explained in Appendix B, in the MC calculations, the 
functional F[z\ is truncated to the following discrete sum 

F'[z] ^ y E «? ^((^ - J>- ^^ - ^i)- (4.11) 

A procedure is introduced then to correct for the missing diagonal {i = j) terms, assuming 
that these are approximately given by the corresponding terms in the calculation of the 
known function F{t). This amounts to correct the raw data by the quantity 

(e^'[-.-])^(e^'[^'-])e^, (4.12) 

with 

C = Hr) - y E «' ^((^ - J>- 0)- (4-13) 

This correction is linear in r, and affects for instance the calculation of the free energy, 
given by the correlator evaluated at r = /3. Table 4.2 summarizes the results. As we see, 
the correction is small and is no more than a few percent. 
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T 


F^^{P,p = 0) 


C 


0.75 


0.382 


0.0102 


1 


0.366 


0.0131 


1.5 


0.346 


0.0184 


2 


0.331 


0.0235 



Table 2: F^'~^{l3,p — 0) for varfous temperatures obtained with the raw MC data and the correction 
C in Eq. (4.12). The values of the free-energy shift obtained here for T=l and 1.5 can be compared 
with the one-loop results in Table 3.1 for the cases T/M = 0.133 and 0.2, respectively. In absolute 
values, the MC (one-loop) free energy shifts are 0.366 (0.336) and 0.346 (0.318), respectively for 
the two cases; the one- loop approximation underestimates the free energy shift. 
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Figure 11: Left panel: -^^^(r, p = 0) obtained with the raw MC data, for all the temperatures 
covered by our analysis. As usual M = 7.5. Right panel: the same quantity after the correction 
indicated in Eq. (4.12) (and labelles here as "MC-t-renorm"). 



The quantity -F (r, p = 0) is shown in Fig. 11 for various temperatures. The left panel 
displays .Fg (t") = -F (r)— (t//3)-F (/3). The curves are obtained employing directly the 
raw MC data, which are not affected by the correction (4.12). In the right panel we show 
the corrected results. As already found in studying the r = correlator, the curves move 
closer to the static result as the ratio T jM decreases, due to the suppression of diffusion. 

In Fig. 12 we compare the MC results with those of the one-loop calculation presented 
in Sect. 3, in which the Euclidean correlator was obtained through the numerical integration 
of the corresponding spectral function, according to Eq. (2.9). The MC points start quite 
close to the one-loop curves corresponding to the same value of T/M, in agreement with the 
expectation that the short-time behavior is governed by perturbation theory (as already 
discussed in Sect. 3 commenting Figs. 8 and 9). For large values of r//3 the MC results 
lie above the one-loop curves. This general behavior is also analyzed within the simple 
toy-model presented in Appendix A. 
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Figure 12: A comparison between -F''^^(r, p = 0) and F^^'''>{t,p — 0) as a function of the ratio 
T/M . The one-loop curves are obtained from a numerical integration of the charm (M = 1.5 GeV) 
spectral density studied in Sect. 3 for T = 200 and T = 300 MeV. In the right panel the set of curves 
corresponding to T = 300 MeV has been translated downwards by -0.005 in order to make the 
figure more readable. 



4.3 The spectral function 

We turn now to the reconstruction of the heavy quark spectral density from the Euclidean 
correlator obtained with the MC calculation. To do so, we need to invert Eq. (2.9), a well 
known difficult problem. We use here a maximum entropy analysis (MEM), according to 
the algorithm described in Ref. [24]. Another exhaustive introduction to the method can 
be found in Ref. [25]. In such an approach, one determines the "best" possible spectral 
function, given the information one has about the Euclidean correlator (the "data"), and 
prior information one has about the spectral density, such as the fact that it is positive 
definite (and hence can be interpreted as a probability density) and that it satisfies some 
sum rules. The procedure involves the maximization of an entropy function (actually the 
minimization of a free-energy), which is defined with respect to a default model: in the 
absence of data, the spectral density coming out of the entropy maximization is the default 
model. There is, of course, a delicate interplay between the effect of the data and that of the 
default model , and the resulting spectral density will in general keep some reminiscence 
of the chosen default model. In order to explore such a systematic uncertainty we will 
consider two different default models: a constant (within the finite range \uj\ < 10), and a 
Gaussian of the form exp[— (cj— M)^/27^]/y27r72. In both cases we adjust the parameters 
of the default model so that the first two sum rules in Eq. (2.16) are fulfilled. 

Throughout this section the results will be expressed in terms of dimensionless vari- 
ables, displaying for instance a{uj)T as a function of a)/T = {uj—M)/T, the only parameter 
left being the ratio T/M. It is then useful to recall that in the static limit T/M — t- 0, the 
free-energy shift is —aniD/'^T ^ —0.476, while Wpi/T ~ 1.373 controls the location of the 
plasmon absorption/emission peaks . 

As the first test of the potentiality of the MEM procedure, and of the systematic uncer- 
tainties attached to the choice of the default model, we reconstruct the (known) one-loop 
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Figure 13: Left panel: a test of the MEM reconstruction of ttie one-loop spectral function for a 
charm quark (M = 1.5 GeV) at three different temperatures. A Gaussian prior is used. The shift 
of the main peak is systematically underestimated, and its width overestimated. Right panel: the 
dependence on the default model. We use a Gaussian and a constant, as explained in the text. The 
MEM procedure was applied to the data for G^'~^{t) in the case T/M = 0.133 and to the exact 
result for G*^^°°(r). The constant default model leads systematically to larger shift and a bigger 
width than the Gaussian default model. 



spectral density from the one-loop Euclidean correlator G^^^' (r) obtained in Sect. 3, through 
the integration of the corresponding spectral density. We use a large set of data points 
(~ 500), and take a heavy quark mass M = 1.5 GeV and temperatures r = 200,300 and 500 
MeV. As one can see on the left panel of Fig. 13, the MEM inversion - here performed 
with a Gaussian prior - is able to identify the main peak. However this is broader and less 
shifted then the exact result: the shift is ~— 0.15, while it is ~— 0.35 in the original one- 
loop spectral density. The method also reconstructs a low-energy secondary bump, though 
less pronounced than the plasmon-absorption peak in the original one-loop spectral func- 
tion, and it appears also at lower frequency (~ —2 compared to —1.45). On the other hand 
no signature of the high-energy secondary peak present in a^^^'{uj) is visible in the MEM 
spectral density. In the right panel of Fig. 13 we illustrate the sensitivity to the choice of 
the default model. There, the MC data at T = 1 are used, as well as the known infinite 
mass correlator G ^^{t). We consider a Gaussian (with various values of the width) and 
a constant prior. With a Gaussian prior with width j/T= 0.1, 0.2 and 0.3, the main peak 
is shifted respectively by, -0.05, -0.15 and -0.15. The presence of a spectral strength at low 
energy seems to be a quite robust feature of the spectrum, though the broader the default 
model, the less pronounced the secondary bump is. In particular for a flat prior one finds 
simply a very large broadening and negative shift of the main peak. 

In Fig. 14 we show the results of the MEM inversion of the MC data for G^*-'(t), for 
various values of T/M, including the exact infinite mass limit corresponding to T/M = 0. 
The left panel corresponds to a Gaussian default model with 7 = aT/2. The resulting 
spectral densities present a broad main peak, slightly shifted with respect to its position in 
the vacuum (M) by an amount roughly proportional to T (the curves in the dimensionless 
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Figure 14: The MEM spectral densities a^'-^{uj) for (different values of T/M. In the left /right 
panel a Gaussian/constant default model is employed. For comparison the curves obtained from 
G^^°°{t) are also shown. The dashed vertical lines correspond to the static free-energy shift 
—amD/2T — 0.476. The dot-dashed vertical lines signal icjpi, where uipi/T = 1.373 is the plasma 
frequency. As clearly seen in the left panel, the Gaussian prior leads to an underestimate of the 
shift of the main peak (here estimated as the static free energy shift), together with an overestimate 
of that of the secondary peak at low energy (here estimated by —ojpi). The dependence on T/M is 
very weak. On the right panel one sees that the dependence on T/M is larger with the constant 
prior, and in line with what one expects (the curves move gradually towards that corresponding to 
the infinite mass limit as T/M decreases). 



units employed lie almost on top of each others), but smaller (by a factor ~ 5) than the 
static free-energy shift —am£)/2T~ 0.476. A secondary low-energy bump, more and more 
displaced with respect to the main peak as T/M decreases, is also visible. In the right panel 
the same data are analyzed using a constant default model. In such a case the spectral 
function exhibits only a broad peak with a sizable negative shift which, as T/M— )-0, results 
~50% larger than the static free-energy shift. Furthermore the MEM spectral density, with 
this choice for the prior, has also a long high-energy tail, at variance with what is found 
with the Gaussian default model. 

Finally in Fig. 15 a comparison between the MEM and the one-loop spectral functions 
is given. The main features discussed above can be seen. In particular, the dependence on 
the default model is striking. A Gaussian leads to a very small shift of the main peak. On 
the other hand the constant default model yields a broader and more shifted peak, whose 
strength extends to low energy till displaying a partial overlap with the secondary bump 
of the curve obtained with the Gaussian prior. 

5. Conclusions 

In this paper, we have presented an approach to the dynamics of heavy quarks in a hot 
plasma based on a path integral for non relativistic particles with a non local (in space and 
time) self-interaction that summarizes the effects of the medium to which the heavy quark 
is coupled. 
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Figure 15: The MEM spectral densities a^'-^{uj) obtained with the two different default models 
(dotted and dot-dashed curves), compared to cr(^')(a;) (continuous curves), for two values oiT/M. 



The path integral providing the heavy-quark Euchdean correlator was evaluated nu- 
merically using Monte Carlo techniques. The results of this numerical evaluation were 
analyzed and compared to those of the one-loop calculation, and to those of an exact 
evaluation of the path integral in the infinite mass limit. We showed that the effect of in- 
teractions is to favor the contribution of straight paths in the path integral, and are indeed 
maximum in the infinite mass limit, where the heavy particle stays at rest. Calculations 
were done for a value of the coupling constant that would correspond in QCD to a value of 
the strong coupling constant a^ ss .3. For such a value the one-loop approximation provides 
a reasonable first approximation, but deviations with the exact Monte Carlo results were 
obtained. The Monte Carlo results move towards those of the infinite mass limit as the 
ratio M/T increases, as expected. 

The Monte Carlo calculations of the Euclidean correlator were performed in coordinate 
space, but a simple integration over the spatial coordinates gave the correlator for zero- 
momentum. This allowed us, in particular, to estimate the shift in the free energy of the 
system that is caused by the addition of the heavy quark. We also used the corresponding 
Euclidean correlator to reconstruct the spectral function, through a MEM analysis. Within 
our implementation of this method, we were able only to reproduce the main qualitative 
features, namely a broad main peak, whose shift is only given semi-quantitatively. A 
secondary structure below the main peak, somewhat reminiscent of the plasmon-absorption 
peak of the one-loop spectral function is also seen, but no secondary structure above the 
main peak is detected, only a long tail at large frequencies is observed (and only with a 
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constant prior). The large sensitivity of the MEM analysis to our choices of default models 
does not allow us to draw more robust quantitative conclusions at this stage. On the other 
hand, the qualitative features that we were able to reconstruct may be enough to draw 
conclusions in the two particle problem, which is our ultimate goal. 

The thorough analysis of the one-particle case that we have presented in this paper 
paves the way for several extensions. Clearly the calculation can be improved in several 
places, and the general setting brought closer to QCD without too much efforts. For in- 
stance, we have seen that the HTL approximation used in the description of the hot plasma 
leads to a somewhat unrealistic description of the effects of collisions. While this affects 
only mildly the heavy quark correlator, and only at small times where the calculation is 
in complete control (being essentially perturbation theory), this feature can be improved 
without changing the basic structure of the problem. The calculation of the Euclidean cor- 
relator of a heavy quark-antiquark pair is within reach. The reconstruction of the spectral 
density of a heavy quark pair from its Euclidean correlator faces the same difficulty as met 
in lattice QCD: on the one hand, this offers opportunities for more detailed comparisons 
between the two approaches, on the other hand we note that our path integral for the 
Euclidean correlator can be calculated with high precision, which could be exploited to de- 
velop new methods of reconstruction of the spectral density. Finally one may contemplate 
the possibility of calculating the path integral directly in real time, perhaps at the cost of 
additional approximations. That would allow us to bypass the problem of the analytical 
continuation, and would open the possibility of numerous applications. 
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A. An exactly solvable toy-model 

In this section we present a toy model illustrating some of the features of the calculations 
that are presented in the main text, in particular those features that emerge in the infinite 
mass limit. The model consists of a fermion of mass M coupled to a single harmonic 
oscillator that represents the "medium" . The Hamiltonian of the system is written as 

where ip'^ and ifj are the creation and the annihilation operators of the fermion, {^, ^^} = 1, 
(j) and vr are respectively the coordinate of the oscillator and its conjugate momentum, 
[(;^, vr] = i, and a^,a the associated creation and annihilation operators, [a, a^] = 1. Since 
[H, if)' if)] = 0, the eigenstates of H can be classified in sectors characterized by the eigen- 
value of the fermion number operator ip'^ij). Since the fermion has no internal degree of 
freedom there are only two sectors to consider, those with V'^V' = a-nd with ip'^ij) = 1. 
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The first sector corresponds to the medium without the fermion, and the Hamiltonian is 
simply that of the oscillator 

Ho = mo (a) a + I/2) . (A.2) 



The sector with il^^il) = 1 mimics the case in which one adds the fermion into the medium. 
The corresponding Hamiltonian reads 

Hi = M + Ho + g^ = M + Ho + V, (A.3) 

and it has the structure of Eq. (2.2). It is easily diagonalized by introducing the shifted 
operators 

b=a+^i= and b^ = a^ + ^^, with [b,b^] = l, (A.4) 



2m^ w2m|, 



so that 



«' = (''"^)+'"''(''"'+0- <"' 



The spectrum of Hi is identical to that of Hq, and the shift in the ground-state energy is 
given by: 



9' _.. _ 9' 



A£ = £i-£o = M - ^f-^ = M -airiD, a = ^^, (A.6) 



2mf) 2mjj 



where we have introduced the dimensionless coupling constant a. The ground state of Hi 
is a coherent state characterized by a non- vanishing expectation value of the field </>: 



(0) = -^ = -V^J—. (A.7) 

mf) V mo 

This expectation value plays the role of the classical field ^0 associated with the polarization 
cloud around the heavy quark. 

One may also consider the non-equilibrium situation that corresponds to adding the 
fermion into the system in its ground state at t = 0. Following this initial perturbation, the 
whole system evolves in time with the Hamiltonian Hi. It is then not difficult to establish 
that the expectation value of (p oscillates around its equilibrium value (A.7) according to 

{(l))t = {(/)) egicosniDt -I), (A.8) 

where (peg is given by Eq. (A.7). 

This result holds unchanged when the oscillator is in thermal equilibrium at tempera- 
ture T, that is, {(p) is not affected by thermal fluctuations. Similarly, because the spectra of 
Hi and Hq are identical, the contributions of thermal fluctuations cancel in the difference 
of free energies of the systems with and without the fermion, with the result that this 
difference remains equal to the shift in the ground state energy given by Eq. (A.6). 

Consider now the Euclidean correlator 

G(r) ^ G>{-zt) ^ {^PiT)^p\0))o, (A.9) 
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where the expectation value (. . . )o = Tr [e ^ ° . . . ] /Zq is taken over states of the medium 
without the fermion. One has: 



G>(-ir) = (e^^^e-^^^t)g = (e^o^^e-^i^V^)o = e"^^ (^H^r^-{H^+v)r 



(A.IO) 



where, in the last expression, one recognizes the evolution operator in the interaction 
representation, so that one can write: 



G>i-iT) 



-Mr I rr 

e { Tr exp 



-g I dT'(t)i{T') 
Jo 



A simple calculation yields then the exact result: 

G>{-iT 

where 



^-Mr ^F(t) 



F(t) 



9 



3 f-T 



dr' / dT"A(T' - t"). 
2 ./o JO 



(A.ll) 



(A.12) 



(A.13) 



Here A(r) is the Euclidean propagator for the field </), satisfying periodic boundary condi- 
tions (A(0) = A(/3)): 



A(r) 



(T(/),(r)(/)/(0)) 



1 



2m 



D L 



e-™Dkl(l + iV) + e'"^l^liV 



with A^ the statistical factor 



N = 



(A.14) 



(A.15) 



At this point let us note that the model depends on several dimensionful parameters: 
the mass M, which simply shifts the overall spectrum, and plays no role in the dynamics; the 
Debye mass rriD which characterizes the response of the system to an external perturbation, 
such as the addition of the fermion; the coupling constant g and the temperature T. We 
shall systematically express the coupling between the fermion and the oscillator in terms 
of the dimensionless coupling a, as in Eq. (A. 6). A look at the propagator (A.14) reveals 
that m~jj appears there as the natural time scale, while the statistical factor depends on 
mo/T. 

It is sometimes convenient to write F{t) as the sum of a term Fi{t) linear in r and a 
term F2{t) that is symmetric around /3/2: 



Fi{t) 
F2{t) = a 



am£)T, 

cosh(m£)(r — 13/2)) — cosh(/3?7i£)/2) 
sinh/3(m/5/2) 



Clearly, 



Fi(3) = Fi(/3) = amD(3 = ^f3A{ioJn = 0), 



(A.16) 



(A.17) 



so that M — (l//3)Fi(/3) = AFq is the difference of free energies of the systems with and 
without the fermion (see Eq. (A. 6)). The last equality in Eq. A.17 emphasizes that F{/3) 
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is entirely given by the zero Matsubara frequency part of the oscillator propagator (A. 14), 
with 

A{iujn) = / dre*'^"^A(r). (A.18) 

Jo 

The function F2{t) vanishes at r = and r = /3, by construction, and has its minimum at 

r = /3/2, with value -F2(/3/2) = — Qtanh(/3r?T,£)/4). The slope at r = is —arriD, so that 

the linear contributions cancel in F = Fi + F2, in accordance with the general result (see 

Eq. (2.13)). This is also obvious from Eq. (A. 13): the small r behavior starts at order r^. 

At quadratic order, we have 

F(r « m^i) c^ ^{4>^) = \amlT\l + 27V), {(j?) = ^(1 + 2N) = A(t = 0). 

(A.19) 
One recovers the general result between the coefficient of r^ and the fluctuation of (f) (see 
Eq. (2.14)). We shall return to the short time behavior of the correlator shortly. 

We now exploit the analyticity of G^ and move to real time. This will allow us in 
particular to get the large time behavior of G^(t). One gets from Eq. (A. 12) 

G>(t) = e-^^*e^^W, (A.20) 

with 



F{t) = ^ f ds [ ds'D{s - s'), D{s - s') = iA{T = is, t = is'). (A.21) 

2 Jo Jo 



,2 f-t rt 

ds 

'0 Jo 

A simple calculation then yields 



F(.t) = g^ / 2 ^ ^ ' A.22 

with D{uj) the Fourier transform of the time-ordered propagator D{t) (see Eq. (2.44)). The 
large time behavior of the correlator follows then from Eq. (2.47): 

Fit > m^^) = —tD{uj = 0) = amot. (A. 23) 

It is entirely determined by the static response of the medium. The comparison with 
Eq. (A. 6) reveals that —F{t)/t is the interaction contribution to the energy shift caused 
by the addition of the fermion (see also Eq. (A. 17)). This is similar to what happens in 
the case of an infinitely massive quark although, in the latter case, a damping term also 
appears next to the free energy shift. No such term appears here because of the discrete 
nature of the spectrum. 

One can also calculate the spectral function. To do so, it is convenient to start with 
the following explicit expression of the propagator (A.20): 

G>(t) = exp [-a{l+2N)] exp [-i{M -amD)t\ x 

X exp [a (iVe*™«* + (l + iV)e-^™o*)] , (A.24) 
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Figure 16: The spectral function (A. 26), for different values of the coupling a. The vertical lines 
refer to the position of the "main peak" at uj = M— amo- The delta functions have been smeared 
to gaussians with width 7 = aT/2, and the mass mo is adjusted as a function of the temperature, 
m_D = TyfoL. As the coupling grows, the individual peaks are smoothed out, leaving a broad, 
structureless, distribution. 



and expand the last exponential in powers of a. One gets 



Q>{t\ = g-a(l+2Af)g-i(Af-amD)t 



00 n " / \ 

^^Y.[) (^^)VP'"«*(l + iV)"~Pe~^ ("-P)"^i5*. (A.25) 



n=0 p=0 

The Fourier transform is then obtained immediately and reads 



00 „ n 

a sr-^ in 



(j{oj) = 2^e~"(i+2iV) y -7 y {Nf{l+NyP6 {u + amo - {n-2p)mD) , 

^-^ n'. ^-^ \n / 

n=0 p=0 ^^^ 

(A.26) 

where we have set cD = w — M. The above spectral density exhibits an infinite number 
of peaks in one-to-one correspondence with the transitions between the eigenstates of Hi. 
The major peak is located at a) = —amu. The expansion of the spectral density to order 
a has peaks centered at tD = —am£) it kmu, with k = 0,1... K. Hence, the larger 
the coupling, the larger the number of peaks giving a sizable contribution to the spectral 
density. 

The spectral density is displayed in Fig. 16 for a wide range of values of the coupling 
a. In order to make contact with the general discussion of a heavy quark in a plasma, we 
choose rriD = \/aT (this implies among other things that the coupling among the plasma 
grows similarly the coupling between the fermion and the plasma particles). Also, for the 
purpose of illustrating the global behavior of the spectral function, we smear the delta 
functions by replacing them by gaussians of variance 7~aT. At small coupling, individual 
peaks are recognized. For large coupling, the smearing that we have introduced erases 
the individual secondary peaks, leaving a broad distribution which spreads over a larger 
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and larger frequency interval as the coupling grows. Note that the main peak, located at 
u] = —arriD is shifted to lower frequency as a grows, but the spectral strength remains 
centered around co ~ M. This behavior may be understood in terms of the sum rules 
satisfied by the spectral function. 

These sum rules are obtained from the derivatives s of G^ (t) at t = 0: 



dt"- ^ ' 



t=0 -'-oo ^'^ 



oo J- 

Gf'aiQ). (A.27) 



These derivatives are easiy calculated and one gets 

— o{uj) = \, / — a;cr(a;) = 0, 

/oo J- /<oo J- 

— £l;V(w) = Qm|,(l + 2iV), / — Q^ a{Gj) = am\). (A. 28) 

-oo 2vr J_^ 27r 

These sum rules explain why the spectral weight remains centered around w = 0, with a 
width increasing with a, while the last sum rule suggest a somewhat larger strength at 
positive (D than at negative uj. Note that the sum rules that are displayed explicitly here 
are at most linear in the coupling a. The first higher order correction, of order a^, enters 
at the level of the a; sum rule. 

Let us now turn to the one-loop approximation for the time-ordered (or retarded) 
propagator. The one-loop self-energy of the fermion is easily obtained: 



Zjfw -|- If]) = arrir-, h 

uo — m£) + 17] uj + m£) + 17] 



(A.29) 



The poles of S for Q = ±771/5 correspond to the energies of the fermion having emitted or 
absorbed a quantum of the oscillator, which represent the leading processes that take place 
at weak coupling. The inverse retarded propagator reads 

G^'^{u) + i7]) = -u) -i7] + T,{Ld + ir]). (A. 30) 

Thus, the propagator has three poles, at values uii solutions of the equation 

CJ^ -u;ml)[l + ail + 2N)]-am% = 0. (A.31) 

The general behavior of the solutions may be easily inferred from the graph displayed in 
Fig. 17. The propagator may then be written as 



^— "* COi — CO 
I 

with the residues given by 

z-^ = \-dY.ldG)\i. (A.33) 

The spectral function takes the form 

o-(w) = 27r ^ Zj (5(w - a)j) . (A.34) 
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Figure 17: Graphical solution of the equation to — S((I>), with both S and (I> expressed in units 
of mjj (the curves plotted are S(cj/q; and ujJol). The self-energy (divided by a) exhibits poles 
at (Ii = ±771/3 . The intersections with the straight line w/a give the locations of the poles of the 
propagator. There is always one pole close to w = 0. At weak coupling this pole has the largest 
residue (the straight line in the figure corresponds to a = 0.36). When the coupling grows the 
other two poles move away as « ±amjy, and their residue saturate the sum rule, leaving very little 
spectral weight on the pole at w « (which asymptotically moves to w = — mD/(l + 2iV)). Note 
that the intersection of S with the vertical axis yields the exact energy shift, S(w = 0) = —amo- 



It can be verified that, in the weak coupling limit, this coincides with the general expression 
(A. 26) expanded to order a. The Euclidean correlator is easily obtained from the spectral 
function, and reads 



G>{-tT) 



~Mt 



E 



Zie 



q-Mt^fW(t) 



(A.35) 



which defines the function F^^^'{t). From the correlator calculated for t = /3, one deduces 
the one-loop free energy shift 



F(^')(/3) = ln 



E 



Zi.e 



(A.36) 



This is to be compared to the exact value F{p) = amo/T: the one-loop calculation 
underestimates the exact result. 

A comparison between the one-loop and the exact result is offered in Fig. 18. As we did 
earlier, we may decompose F^^^\t) = F^^^'V) + H^^\r), with F^ 'V) = (.r/P)F^^^\f3). 
The function F2 (r) is plotted in the right panel of Fig. 18. The agreement of the exact 
and one-loop correlators may be understood from the fact that the sum rules (A. 28) are 
exactly satisfied at one loop, namely 



^Zi = l, '^ZiUJi = 0, ^ 



ZiUJ^ 



am£, 



(l + 2iV), Y. 



ZiOJi 



am£,. 



(A.37) 



To these we should add the relation Yli'^i — 0' ^^^^ derives immediately from Eq. (A. 31). 
The sum rules (A.37) are the one-loop transcription of the exact sum rules mentioned 
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Figure 18: Left: The function F^^^^{t) (continuous curves) compared to the exact result F{t) 
(dot-dashed curves) as a function of r//3 for large values of the coupling constant a = 1 and a = 2 . 
Even in this strong coupling regime, the one-loop approximation gives an excellent approximation 
to the exact result up to values r ^ (3 /2a. Right: The same for F2 (r) and F2{t). The slope at 
the origin is a measure of the free energy shift, which can be also read on the left panel as the value 
of F{P). Both plots exhibits clearly that the free-energy shift is underestimated in the one- loop 
approximation. Note also the asymmetry, growing with increasing coupling, of the one-loop results 
with respect to r = /3/2, in contrast to the exact curves. 



above, Eq. (A. 28). They hold exactly at one-loop because the small time behavior of the 
propagator involves also a small g expansion and, as it has already been mentioned after 
Eq. (A. 28), up to order r^, the small r expansion involves terms of the weak coupling 
expansion only up to order g^. Such terms are taken into account exactly by the one- loop 
self energy. The fact that the one-loop result is sufficient to describe the small-r behavior 
appears clearly in the left panel of Fig. 18 where, even for large values of the coupling, the 
one-loop and exact correlators are hardly distinguishable for r//3 small enough. 

Finally in Fig. 19 we provide a comparison between the one-loop spectral function, 
given by Eq. (A. 34), and the exact result. For weak coupling they look quite similar. On 
the other hand for larger coupling more and more secondary peaks contribute to the exact 
spectrum, while the one-loop result can display only three peaks. These, having to fulfill 
the sum-rules (A. 37), result largely distorted. 

B. Details on the path integral implementation 

The path integral that we want to evaluate has the following form (see Eq. (4.5)): 

/;pzF[z,T]e-^"[^'^ 



{F[z,r] 



/;^- 



(B.l) 



For any chosen value of r = N-rttr, where A',- is an integer and Ur a fixed time interval, the 
paths are defined by a discrete set of points {z(rj)}, where tj (0 < tj < r) is a multiple of 
Ur- We choose natural units h = c = ks = ^ and fix the unit of length to be 1 fm, and 
correspondingly the unit of energy (or temperature) to be 197.3 MeV. The time step is 
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cc=0.36 (one-loop) 
a-1. 00 (one-loop) 
a^2.00 (one-loop) 
cc-0.36 (exact) 
a=l. 00 (exact) 
a=2.00 (exact) 
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Figure 19: The one-loop spectral function compared to the exact one. The curves refer to different 
values of the coupling a, from weak (left panel) to strong (right panel). The one- loop result 
(continuous curves), characterized by the presence of only three peaks, is compared to the exact 
one (dot-dashed curves), which has a richer structure. In plotting the curves the delta functions 
have been smeared to gaussians in the same way as in Fig. 16. 



fixed at the value a,- = 0.01 fm/c, and the heavy quark mass at M = 7.5 (corresponding 
approximately to the mass of a charm quark). 

The path integral Monte Carlo method is based on the generation of a Markov chain 
that samples a set of paths according to the distribution 



Wa[z] 



exp(-5a[^]) 
IP 



z\ e" 



(B.2) 



Then the average {F[z])a is evaluated as the arithmetic average over the paths generated 
by the equilibrated Markov chain. 

We used the Metropolis algorithm to generate the Markov chain by a sequence of 
elementary moves. A move is defined as follows: starting from a given path z, one selects 
at random a time tj, and displace the corresponding point z^n) by a quantity 5z uniformly 
distributed in a cube of side d centered at z(tj), thus defining a new path z'; the move is 
accepted with the probability 



vr 



mm 



eM-Sa[z']) 
' exp(-5Q,[2;]) 



(B.3) 



In the present calculation we start from a straight path connecting the origin (0, 0) to the 
point {T,r), and perform at least 10^ x Nr moves to reach equilibrium. During this stage 
the value of d is adjusted to keep the acceptance ratio of attempted moves between 0.45 
and 0.55. Typical values of d were found in the range 0.07 — 0.08, at the temperature 
T = 1.0. Once the Markov chain has reached equilibrium, one continues generating paths, 
and the corresponding paths are used in calculating the average values of {F[z])a- At 
least 10^ X Nt paths of the equilibrated chain are used in the calculation of the average 
value. Finally, the integrand appearing in the right hand side of Eq. (4.6) is evaluated on 
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an equispaced array of 10 points in the interval from a = to a = 1. The resulting curve 
is then interpolated with a cubic spline and integrated using the adaptive Gauss-Kronrod 
method as implemented in the GNU Scientific Library [26]. 

At a given temperature we take tipically between 10 and 20 values of r to determine 
G{T,r), with r varying between and 2 fm. Because a large number of paths are used in 
the calculation of the average, the statistical errors are negligible: one gets relative errors 
of order 10~^ for small r and 10"*^ for the largest values of r. Furthermore, since the 
average is taken over a different set of trajectories at each r, the results at various r are 
uncorrelated. 

There is one issue in this calculation that deserves further comments. It concerns the 
calculation of the integral (4.3). The simplest discretized form of this integral reads 

F[z, Nr] ^ y E «' ^((^ - J>- ^^ - ^^■)- (^-4) 

This, however, cannot be used as it stands since, as we have seen in Sect. 2.2, A(0, 0) is 
logarithmically divergent, so that the diagonal terms i = j in the expression above are ill 
defined. Before we explain how we have gone around this difficulty, let us examine the 
calculation of the same integral in the infinite mass limit, where the paths are frozen at 
the origin, i.e., z{Ti) = for all i. Then the functional F[z,r] reduces to the function 
F{t) (Eq. (2.52)), that we may write, using the same discretization as above but for the 
"diagonal" terms, as 

F{t) ~ y E a? A((i - j)ar, 0) + iV,F(a,), (B.5) 

where Fla-j-) is the exact value of the integral on a square of side a,-. The comparison of the 
value of F{T) — NF{ar) calculated exactly, and from the discretized sum in Eq. (B.5), yields 
an estimate of the discretization error in the evaluation of the integral. As can be seen in 
Fig. 20, this error is of (relative) order 10^^ and increases slightly towards small values of 
r. We note also that for t//3 ^ 0.2, the contribution of the diagonal terms amounts to less 
than 10%. We have exploited these features in order to make the following simplifications: 

i) Only the off-diagonal terms are used in the sampling of paths, that is F[z,t] is 
replaced for that purpose by 

ii) A correction is applied to compensate for the omission of the diagonal terms, as- 
suming this correction to be given by NT-F{ar) (in practice we calculate this correction 
from the difference between F{t) and F'[z,r] in Eq. (B.6) estimated for z = 0). Note that 
this correction is presumably an overestimates. Indeed because of diffusion, at time ~ a,-, 
the heavy quark is on the average at a distance f = ^Zt jIM away from the origin, and 
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Figure 20: Left panel: the relative error in evaluating F{t) with the discretization algorithm 
employed for the path integral {Njs = 100). The black dashed curve is obtained by simply dropping 
the i = j terms in the sum, as in Eq. (B.6); the red dot-dashed curve arises after correcting with 
the term NrFiar) from Eq. (B.5). Right panel: The short-timc/distance behavior of the HTL 
propagator A(t, r), for various values of t/(3 (blow-up version of Fig. 1). The vertical lines refer to 
the corresponding values of r= Vv^^T)- 



A(r, f) < A(r, 0). A quantitative measure of this overestimate (which is of the order 30%) 
can be read off the right panel of Fig. 20. 

A final source of errors comes from the fact that the MC calculation has been set in 
fixed physical units. Thus, as we change the temperature, of equivalently /3 = 1/T, one 
varies the number of discretization points, with T = 1 (~ 200 MeV) corresponding to 
Nj3 = 100. Increasing T, means decreasing /3, and correspondingly A^^. Because of this, 
the calculations become less precise as the temperature increases. 
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